Problema 1

Dada una matriz \(4 \times 4\) de variables aleatorias Bernoulli denotada por \(\displaystyle \left[X_{ij}\right]\), sean \(N(X)\) el número total de éxitos en \(X\) y \(D(X)\) el total de vecinos (horizontales o verticales) de dichos éxitos que difieren. Nótese que, si se escriben los elementos de \(X\) en un vector \(V\), existe una matriz \(\displaystyle M_{24 \times 16}\) tal que \(D(X)\) es la suma de valores absolutos de \(MV\). El \(24\) surge porque hay \(24\) pares de vecinos a revisar.

Asimismo, supongamos que la distribución de \(X\), \(\pi(X)\), es proporcional a

\[ \pi(X) \propto p^{N(X)} (1 - p)^{16 - N(X)} \exp \left( -\lambda D(X) \right). \]

Si \(\lambda = 0\), las variables son independientes e idénticamente distribuidas \(\text{Bernoulli}(p)\).

Hay \(2^{16}\) posibles estados, uno por cada valor posible de \(X\). Usaremos el método Metropolis-Hastings con los siguientes kérneles de transición:

  1. Sea \(q_1\) tal que cada transición es igualmente plausible con probabilidad \(1/2^{16}\). Es decir, el siguiente estado candidato para \(X\) es un vector de \(16\) \(\text{Bernoulli}(p)\) iid.

  2. Sea \(q_2\) tal que se elige una de las 16 entradas en \(X\) con probabilidad \(1/16\), y luego se determina el valor de la celda a ser \(0\) o \(1\) con probabilidad \(0.5\) en cada caso. Entonces, a lo más puede cambiar un solo elemento de \(X\) en cada transición.

Ambas \(q_1\) y \(q_2\) son simétricas, irreducibles y tienen diagonales positivas. Es fácil ver que la primera se mueve más rápido que la segunda.

Estamos interesados en la probabilidad de que todos los elementos de la diagonal sean \(1\). Al usar \(q_1\) y \(q_2\), estimaremos dicha probabilidad para los valores \(\lambda \in \{0, 1, 3\}\) y \(p \in \{0.5, 0.8\}\).

Problema 2

Queremos investigar el desempeño del algoritmo de Metropolis-Hastings usando una caminata aleatoria cuando la distribución objetivo es una mezcla de dos densidades normales bivariadas para \(\theta = (\theta_1, \theta_2)\)

\[ \pi(\theta) = 0.7 \mathscr{N} (\theta \mid \mu_1,\ \Sigma_1) + 0.3 \mathscr{N} (\theta \mid \mu_2,\ \Sigma_2) \]

donde

\[ \mu_1 = \begin{pmatrix} 4 \\ 5 \end{pmatrix},\ \mu_2 = \begin{pmatrix} 0.7 \\ 3.5 \end{pmatrix},\ \Sigma_1 = \begin{pmatrix} 1 & 0.7 \\ 0.7 & 1 \end{pmatrix},\ \Sigma_2 = \begin{pmatrix} 1 & -0.7 \\ -0.7 & 1 \end{pmatrix}. \]

Supongamos que no podemos muestrear \(\pi\) directamente y, en su lugar, usamos Metropolis-Hastings con

\[ q(y \mid x) = \mathscr{N}(y \mid x, v I_2) \]

donde \(v\) es un parámetro de ajuste.

mu.1 <- c(4, 5)
mu.2 <- c(0.7, 3.5)
sigma.1 <- matrix(c(1, 0.7, 0.7, 1), 2, 2)
sigma.2 <- matrix(c(1, -0.7, -0.7, 1), 2, 2)
I.2 <- diag(2)

mix.norm <- function(x) {
  0.7 * dmvnorm(x, mean = mu.1, sigma = sigma.1) +
    0.3 * dmvnorm(x, mean = mu.2, sigma = sigma.2)
}

q <- function(a, b, c) {
  dmvnorm(a, mean = b, sigma = (c * I.2))
}

# Simulación
simulation <- function(n, init, v) {
  X <- matrix(c(NA, NA), 2, n)
  X[, 1] <- init
  count <- 0
  
  for (t in 1:(n - 1)) {
    x <- X[, t]
    y <- x + (2 * runif(2) - 1)
    U <- runif(1)
    MHR <- (mix.norm(y) * q(x, y, v)) /
      (mix.norm(x) * q(y, x, v))
    alpha <- min(1, MHR)
    if (U <= alpha) {
      X[, t + 1] <- y
    }
    else {
      X[, t + 1] <- x
    }
    ifelse(alpha != 1,
           count <- count + 1,
           count <- count)
  }
  
  prob <- count / n
  X <- data.frame(x = X[1, ], y = X[2, ], n = 1:n, prob = prob, v = v)
  
  return(X)
}

# Visualización
visualize <- function(X) {
  colours <-
    c("Punto inicial" = "#1e40ca",
      "Media posterior" = "#00a2ed")
  
  p.1 <- ggplot(X) +
    geom_path(aes(x = x, y = y), size = 0.1, alpha = 0.3) +
    geom_point(aes(x = x, y = y), size = 0.1) +
    geom_point(aes(x = x[1],
                   y = y[1],
                   colour = "Punto inicial"), size = 2) +
    geom_point(aes(
      x = mean(x),
      y = mean(y),
      colour = "Media posterior"
    ), size = 2) +
    labs(
      title = NULL,
      x = expression(theta[1]),
      y = expression(theta[2]),
      caption = paste("Tasa de aceptación:", X$prob[1],
                      "\nv =", X$v[1])
    ) +
    scale_colour_manual(name = NULL, values = colours) +
    theme(legend.position = "top")
  
  p.2 <- ggplot(X) +
    geom_line(aes(x = n, y = x)) +
    labs(title = NULL,
         x = "Iteración",
         y = expression(theta[1]))
  
  p.3 <- ggplot(X) +
    geom_line(aes(x = n, y = y)) +
    labs(title = NULL,
         x = "Iteración",
         y = expression(theta[2]))
  
  p.4 <- ggplot(X, aes(x = x, y = ..density..)) +
    geom_histogram(binwidth = 0.5,
                   fill = "#1e40ca",
                   alpha = 0.6) +
    geom_density() +
    labs(title = NULL,
         x = expression(theta[1]))
  
  p.5 <- ggplot(X, aes(x = y, y = ..density..)) +
    geom_histogram(binwidth = 0.5,
                   fill = "#1e40ca",
                   alpha = 0.6) +
    geom_density() +
    labs(title = NULL,
         x = expression(theta[2]))
  
  plots <- list(p.1, p.2, p.3, p.4, p.5)
  return(plots)
}

Al momento de otorgar valores pequeños o grandes a \(v\) en \(q(y \mid x)\), el efecto de dicha \(v\) es nulo pues, al crear nuestra

\[ \alpha (x, y) = \min \left\{1, \frac{\pi(y) q(x \mid y)}{\pi(x) q(y \mid x)}\right\}, \]

tanto el numerador como el denominador relacionados a nuestra razón de Hastings se escalan a la misma varianza.

Nótese lo anteriormente mencionado en las siguientes caminatas aleatorias para \(v = 0.01\) y \(v = 100\):

set.seed(181121)
visualize(simulation(1000, c(0, 0), 0.01))[[1]]

set.seed(181121)
visualize(simulation(1000, c(0, 0), 100))[[1]]

Asimismo, podemos visualizar la función de autocorrelación de \(\theta_2\):

set.seed(181121)
Y <- simulation(1000, c(0, 0), 1)
acf(Y$y, main = expression(theta[2]))

Finalmente, con \(5,000\) extracciones, nótese que nuestra tasa de aceptación es de \(0.6802\). Gráficamente:

set.seed(181121)
Y <- simulation(5000, c(0, 0), 1)
visualize(Y)[[1]]

visualize(Y)[[2]]

visualize(Y)[[3]]

visualize(Y)[[4]]

visualize(Y)[[5]]

Problema 3

Sea \(\displaystyle \theta = \int_0^{\pi/3} \sin(t)dt\), usaremos Monte Carlo para calcular un estimador \(\hat{\theta}\).

Monte Carlo crudo

set.seed(181121)
n <- 100000
a <- 0
b <- pi / 3

sim.u     <- runif(n, a, b)
sim.raw   <- sin(sim.u)
theta.raw <- cumsum(sim.raw) / (1:n)
theta.raw <- (b - a) * theta.raw
theta.1   <- tail(theta.raw, 1)
var.1     <- var(sim.raw)

Nuestro estimado es: \(\hat{\theta}_1 = 0.5003\).

Variables antitéticas

sim.v     <- a + (b - sim.u)
sim.ant   <- (sin(sim.u) + sin(sim.v)) / 2
theta.ant <- cumsum(sim.ant) / (1:n)
theta.ant <- (b - a) * theta.ant
theta.2   <- tail(theta.ant, 1)
var.2     <- var(sim.ant)

Nuestro estimado es: \(\hat{\theta}_2 = 0.4999\). Asimismo, con variables antitéticas, la varianza se reduce en un \(99.3869\%\).

Variables de control

Sabemos \(\sin(t) = \cos(t - \pi / 2)\). Definamos \(y = t - \pi / 2\) como variable de control.

k <- 1000

# Piloto
set.seed(181121)
u <- runif(k, a, b)
x <- pi / 2 - u
y <- cos(x)
alpha <- -lm(y ~ x)$coeff[2]

# Simulación
set.seed(181121)
u <- runif(n, a, b)
x <- pi / 2 - u
y <- cos(x)

sim.cv   <- y + alpha * (x - mean(x))
theta.cv <- cumsum(sim.cv) / (1:n)
theta.cv <- (b - a) * theta.cv
theta.3  <- tail(theta.cv, 1)
var.3    <- var(sim.cv)

Nuestro estimado es: \(\hat{\theta}_3 = 0.5003\). De tal forma que, con nuestra variable de control, la varianza se reduce en un \(99.3659\%\). Finalmente, nuestro estimador \(\hat{\theta}_3\) reduce su varianza en un \(3.4153\%\) respecto a \(\hat{\theta}_2\).

Gráficamente,

Problema 4

Dados los siguientes números normales con \(\mu = 3\) y \(\sigma = 1\)

obs <- c(4.59, 4.153, 2.46, 2.732, 2.973)

la función de distribución empírica \(F_n(x)\) está dada por:

\[ F_5(x) = \begin{cases} 0, & \text{si}\ \ x_{(1)} > x \\ 1/5, & \text{si}\ \ x_{(1)} \leq x < x_{(2)} \\ 2/5, & \text{si}\ \ x_{(2)} \leq x < x_{(3)} \\ 3/5, & \text{si}\ \ x_{(3)} \leq x < x_{(4)} \\ 4/5, & \text{si}\ \ x_{(4)} \leq x < x_{(5)} \\ 1, & \text{si}\ \ x_{(5)} \leq x \end{cases} \]

Asimismo,

\[ \begin{align} P(F_5(3) \leq 2/5) &= P(F_5(3) = 0) + P(F_5(3) = 1/5) + P(F_5(3) = 2/5) \\ &= {5 \choose 0} F(3)^0 \left( 1 - F(3) \right)^{5 - 0} \\ &+ {5 \choose 1} F(3)^1 \left( 1 - F(3) \right)^{5 - 1} \\ &+ {5 \choose 2} F(3)^2 \left( 1 - F(3) \right)^{5 - 2} \end{align} \]

donde \(F(x)\) es una normal con \(\mu = 3\) y \(\sigma = 1\).

Así pues,

mean <- 3
sd <- 1

prob <- 0

for (i in 0:2) {
  p <- choose(5, i) *
    (pnorm(3, mean, sd)) ^ (i) *
    (1 - pnorm(3, mean, sd)) ^ (5 - i)
  prob <- prob + p
}

prob
## [1] 0.5

De forma analítica, es fácil ver que \(F(3) = \Phi(0) = 0.5\) donde \(\Phi(\cdot)\) es una normal estándar. De igual modo, nótese que \(F_5(x) \leq 2/5\) si y solo si \(x < x_{(3)}\) donde \(x_{(3)}\) es la medida de nuestra muestra. Por lo tanto, podemos concluir que \(P(F_5(3) \leq 2/5) = 0.5\).